function [assignments] = FMC_interpret(AllSals, numClusters, AllPs, A_D, beta)
%% Interpretation of FMC data and grain assignment

%% run interpret 3D


A_D = A_D+A_D';
NNlist = {};
for i = 1:length(A_D)
    NNlist{i} = find(A_D(i,:));
end

%% sort by saliency
% sort all course nodes for all s by saliency from lowest to highest
[~,elements] = sort(AllSals,'ascend'); %elements is list of indices representing the sorted AllSals list
% all points assigned to own node Nx2 [assignment,confidence]
% pointData is euler1 euler2 euler3 x y z confidence
numPoints = numClusters(1);
assignmentsX=zeros(numPoints,1);  % assignment of each poing
assignmentsY=zeros(numPoints,1);  % probability of that assignment

%% s values of each node
% vector sAssign has s value of all nodes
numS = length(numClusters); %number of iterations of s
sAssign=ones(numClusters(1),1); %initiating sAssign
lengthsub=zeros(1,numS);


%1 spot in sAssign for each node, increasing number for each saliency set
parfor curS=2:numS
    sAssign=[sAssign; curS*ones(numClusters(curS),1)];    
end

% build assignments
for curS = numS-4:numS
    fprintf('Considering S = %i \n', curS)
    
    %find nodes for current s
    elemS = elements(sAssign(elements) == curS)';   %nodes of curS
    sLen = numClusters(curS);
    numNodes = length(elemS);
    allUs = zeros(numPoints,numNodes+1);
    
  
    % for each s from current s to s=1, iteratively find which nodes
    % are inside of the node by multiplying by Ps
    backPs = AllPs{curS};
    for iterS=curS-1:-1:2
        backPs=AllPs{iterS}*backPs;
    end
    
    % find probability of points assigned to each node
    parfor curNodeIndex = 1:numNodes
        %U is a vector of zeros whose entries correspond to whether
        %subnodes belong to the current node...
        U=zeros(sLen,1);
        U(curNodeIndex)=1;  %seeded with a one for the current node
        U = backPs*U;
        
       allUs(:,curNodeIndex) = U; 
        
    end
    
    allUs(:,numNodes + 1) = assignmentsY;
       
    
    %assign clusters to points based on the strongest probability of any
    %point in the cluster
    for point = 1:numPoints
        [Prob, NodeIndex] = max(allUs(point,:));
        if ((NodeIndex <= numNodes) && (Prob ~= 1))
            nodePoints = allUs(:,NodeIndex)>=beta;
            assignmentsX(nodePoints) = NodeIndex;
            assignmentsY(nodePoints) = Prob;
        end
    end
end


%% Break up non-continuous (Michigan) clusters
rLimit = 15;   %Recursion limit (MATLAB tends to crash around 840)
%N.B. This runs considerably faster with a low recursion limit due to
%MATLAB's huge overhead. Better to make more calls than fewer, deeper calls
% set(0,'RecursionLimit',rLimit)

disp('Breaking non-continuous clusters')

flag = zeros(length(assignmentsX),1);
assignmentsN = zeros(length(assignmentsX),1);
contPoints = [];
for point = 2:length(flag)
    if assignmentsN(point) == 0
        oldClust = assignmentsX(point);
        newClust = point;
        if (flag(point)>1)     %if the point has been reserved, continue the same grouping
            newClust = flag(point);
        end
        assignmentsN(point) = newClust;
        rDepth = 1;

        [assignmentsN, flag, contPoints] = clusterBreaker(assignmentsN, point, assignmentsX, flag, NNlist, oldClust, newClust, contPoints, rDepth, rLimit-5);
        
        while isempty(contPoints) == 0
            pointC = contPoints(1);
            contPoints = contPoints(2:end);
            assignmentsN(pointC) = newClust;
            rDepth = 1;
            [assignmentsN, flag, contPoints] = clusterBreaker(assignmentsN, pointC, assignmentsX, flag, NNlist, oldClust, newClust, contPoints, rDepth, rLimit-5);
        end
    end
end

singlePoints = find(flag == -1);
for pt = singlePoints'
assignmentsN(pt) = assignmentsN(NNlist{pt}(2));
end

assignmentsX = assignmentsN;
set(0,'RecursionLimit',500)


assignments = [assignmentsX assignmentsY];
    
end
